%% mu
%unit: Pa*s
function mu=mu_cal(T,P)
load mu_data
T=T+273.15;
[~,c]=min(abs(mu_T-T));
if P>mu_P(end)
    P=mu_P(end);
elseif P<mu_P(1)
    P=mu_P(1);
end

mu=interp1(mu_P',mu_data(:,c),P);
end